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1. Introduction 

An evolutionary (time-dependent) partial differential equation (PDE) can be reduced to a system 
of ordinary differential equations (ODEs) by replacing the spatial derivatives with finite-difference 
approximations. The resulting approximation is called semidiscrete since the time variable is left 
continuous. The procedure of reducing a PDE to an ODE system is often called the method of lines 
since a solution of the ODE system gives an approximation to the PDE solution along x equals 
constant lines in (x, t) space. 

A semidiscrete approximation for a hyperbolic initial-boundary-value problem (IBVP) leads 
to a complication since, in general, more boundary conditions are required for the semidiscrete 
approximation than are specified for the PDE. The additional boundary conditions are often called 
numerical boundary conditions. Any numerical procedure used to provide a numerical boundary 
condition is called a numerical boundary scheme (NBS). 

An essential requirement for a semidiscrete approximation is the convergence (as the spatial 
mesh is refined) of the approximate solution to the solution of the PDE. If an approximation is 
consistent, then by the Lax equivalence theorem stability is a necessary and sufficient condition 
for convergence. It is generally easy to check the consistency of the approximation; however, the 
stability analysis can be a formidable problem. 

Improper treatment of the NBS can lead to instability of the semidiscrete approximation even 
though one starts with a stable approximation for the pure initial-value problem (IVP) or Cauchy 
problem. For the purposes of this paper, we will assume that the approximation is consistent and 
stable (i.e. convergent) for the IVP and will consider only the effect of the NBS on the stability of 
the semidiscrete IBVP. 

For a linear IVP or an IBVP with homogeneous boundary conditions, a semidiscrete approxi- 
mation (on J spatial mesh intervals) results in a system of ODEs of the form du(t)/dt = Au(<) 
where u is a J-component vector, and A is a J x J matrix. The solution of the homogeneous 
ODE system can be written as u(t) = e At u(0) where e At is an exponential matrix. The classical 
Lax- Rich trnyer stability definition is equivalent to the requirement that the matrix norm of e At 
be uniformly bounded for 0 < t < T independent of the spatial mesh size. The practical problem 
is that there is no known simple algebraic test for the uniform boundedness of the matrix norm 
for hyperbolic IBVPs. Necessary (and sufficient in special cases, e.g., IVPs) conditions can be ex- 
pressed in terms of the eigenvalues of the matrix A; however, for the IBVP the eigenvalue analysis 
is, in general, intractable and the resulting conditions are not sufficient for stability. 

In the 1960s and early 1970s, a stability theory for fully discrete IBVPs was developed by 
Godunov and Ryabenkii [l], Kreiss [3], Osher [4], and Gustafsson, Kreiss, and Sundstrom [2]. For 
the purposes of this paper, we refer to this theory as the GKS theory. Trefethen [6] showed that the 
main result of the theory has a physical interpretation in terms of group velocity, and Strikwerda 
[5] extended the theory to semidiscrete approximations. 

The GKS theory reduces the stability analysis of a semidiscrete IBVP on a finite domain to the 
study of three auxiliary problems: the Cauchy problem and the related left- and right-quarter- 
plane problems. The stability of the left- and right-quarter-plane problems is checked by the 
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normal mode analysis. An advantage of the GKS theory is that it provides an algebraic test which 
is necessary and sufficient for stability but does so at the expense of using a more comnlicated 
norm (in the stability definition). 

Since the Lax-Richtmyer and GKS stability definitions differ, the connection between the eigen- 
value analysis for the finite domain problem (which gives necessary Lax-Richtmyer stability con- 
ditions) and the eigenvalue or normal mode analysis for the quarter plane problems (which gives 
necessary and sufficient GKS stability conditions) is rather obscure. In this paper, we consider a 
direct algebraic comparison between the stability polynomial of the finite domain problem and the 
polynomials associated with the quarter-plane problems of the GKS theory. We show under what 
asymptotic conditions the finite-domain analysis leads to a stability polynomial which is consis- 
tent with the quarter plane analysis. The asymptotic eigenvalue analysis (finite-domain, J — ► oo) 
establishes the connection between the algebraic tests of the GKS stability theory and the Lax- 
Richtmyer stability definition. In addition, it also leads to a conjecture which gives necessary and 
sufficient conditions for Lax-Richtmyer stability in terms of the algebraic tests of the GKS theory. 

This short paper gives only a brief treatment of our analysis. A detailed exposition (including 
plots of eigenvalue distributions for various NBSs) is given by the authors in [8]. This paper is 
restricted to semidiscrete approximations. For fully discrete difference approximations we have 
recently stated a conjecture [7] which relates the Lax-Richtmyer and GKS stability theories by a 
generalization of a theorem from linear algebra. 

2. Initial-Boundary- Value Problem for a Model Hyperbolic Equation 

For simplicity we restrict our attention to the stability of semidiscrete approximations to the 
IBVP for the model hyperbolic equation 

du du . , 

^=c— , 0 < x < L, t>0 (2.1) 

where c is a real constant. Initial data u(x,0) = /(x), 0 < x < L are given at t = 0 and for c > 0 
the problem is well-posed if an analytical boundary condition is prescribed at x — L 

u(L,t) = y(t) for c > 0. (2.2) 

3. A Prototype Semidiscrete Approximation for the Model IBVP 

Let the spatial interval L be divided into J subintervals of length Ax, i.e., J Ax = L, x = xy = 
j Ax. On the interior mesh points the spatial derivative u x is replaced by a second order central 
difference quotient and we obtain the system of ODEs 

^■ = 2 ^ (Uy+1_Uy - l) ’ i = 1 ’ 2 ’ " >J_1 (31) 

where uy(t) = uy denotes the semidiscrete approximation to u(x,t). The right boundary (x = L) 
is advanced by using the analytical boundary condition (2.2). We assume that the boundary 
condition is homogeneous, i.e., g(t) = 0 and for the semidiscrete problem we write 


uj = 0. (3.2) 

The spatial computational stencil of (3.1) uses the 3 points j,j ± 1. If we apply (3.1) at the left 
boundary (j = 0), then the stencil protrudes one point to the left of the boundary. It is clear 
that an additional numerical boundary scheme (NBS) is required to determine the semidiscrete 
solution. At this boundary [j = 0) we change from a centered approximation to a one-sided 
spatial differencing approximation of u x : 

= ~~[-au 2 + (1 + 2a)ui - (1 + a)u 0 ] (3.3) 

o ** x 
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du 

dx 


where a is a parameter. The approximation (3.3) is first-order accurate for any a except a = 1/2 
in which case it is second-order accurate. If we insert (3.3) into the PDE (2.1) evaluated at j = 0, 
there follows the NBS 

duo c 

= ^[ -aU2 + (* + 2a ) u i ~ (1 + a)«o)- (3.4) 

The system of ODEs (3.1) together with the analytical boundary condition (3.2) and the NBS 
(3.4) can be written in vector-matrix form as 


du(f) 

It 


— Au(f) 


(3.5) 


where u is a ./-component vector and A is a J x J matrix. 

4 . Lax-Richtmyer Stability of a Semidiscrete Approximation 

The essential element in the stability of a semidiscrete approximation represented by (3.5) is the 
behavior of the solution as the spatial mesh is refined (Ax — ► 0, or J — ► oo). Consequently, one 
must consider an infinite sequence of ODE systems. The J-th member of the sequence is the ODE 
system (3.5) of dimension J. In order to define stability, we need some measure of the magnitude 
of the solution vector and we use a conventional vector norm || • || . 

Lax-Richtmyer stability for a semidiscrete approximation is defined as follows: 

Definition 4 . 1 . A semidiscrete approximation represented by the sequence of ODEs (3.5) is said 
to be Lax-Richtmyer stable if there exists a constant K > 0 such that for any initial condition u(0) 


M*)0<*Ho)|| 


(4.1) 


for all J, J Ax = L with L fixed and for all t , 0 < t < T with T fixed. 

The eigenvalues st of A are, in general, complex and we write st = 9i(s/) + t‘9f(s*). Any 
semidiscrete approximation for (2.1) will have the factor c/Ai on the right-hand side (see e.g. 
(3.1)). In considering the eigenvalues of the matrix A it is convenient to define a matrix A such 
that A = ( | c | / Ax) A. Consequently, the eigenvalues I of A are related to the eigenvalues s of A by 
s = sAx/|c|. 

A necessary condition for Lax-Richtmyer stability can be stated as follows: 

Lemma 4.1. A necessary condition for a semidiscrete approximation to be Lax-Richtmyer stable 
is that for all J, J Ax = L and all t, 0 < t < T there exists a nonnegative constant w such that 

maxUt(se) < y . (4.2) 


If the matrix A is a normal matrix, i.e., AA* = A* A, then (4.2) is both necessary and sufficient 
for stability in the L 2 norm. It is important to realize that a stable semidiscrete approximation 
can have eigenvalues st with positive real parts, but the real parts must approach zero at least as 
fast as 1/J. 

5. Normal Mode Analysis (Quarter-Plane Problems) 

The GKS theory (see introduction) reduces the stability analysis of an IBVP on a finite domain 
to the study of the three auxiliary problems: the Cauchy problem and the related left- and right- 
quarter-plane problems. For example, one can obtain the related right-quarter-plane problem from 
the finite domain semidiscrete problem by fixing Ax and the left boundary at x = 0 and letting 
J —* 00 . Note that now L is not fixed but L — ► 00 as J — *■ 00 , and the resulting spatial domain is 
[0 < x < 00 ). 

The algebraic tests of the GKS theory are carried out by means of the normal mode analysis 
which is based on the resolvent equations. The resolvent equations for a quarter-plane problem 
are obtained by substituting 

uj(t) = e‘%-, (5.1) 
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where s is a complex constant, into the semidiscrete approximation (i.e., the interior scheme) and 
the boundary conditions. The resolvent equations consist of a difference equation and boundary 
conditions for the eigenfunction <f>j. The general solution of the resolvent equations which is in Li 
has the form 

<f>j = <£ 0 /c J , |/c| < 1. (5.2) 

The semidiscrete approximation is GKS stable if there are no nontrivial solutions of the form 


5R (s) > 0, \k\ < 1 


(5.3) 


and 

5R(s) =0, |ac| = 1 such that if |/c*| — ► 1“ then 5R(s*) —* 0 + (5.4) 

where k* indicates a perturbation off the unit circle and s* a perturbation off the imaginary axis 
(complex s-plane). We refer to an eigenvalue of the form (5.3) or (5.4) as a GKS eigenvalue or a 
GKS generalized eigenvalue, respectively. 

To illustrate the application of the quarter-plane normal mode analysis we compute the right 
quarter-plane eigensolutions for the semidiscrete approximation described in section 3. If we sub- 
stitute (5.1) into the interior scheme (3.1) and into the NBS (3.4) we obtain 

2sAx sAx 

<f>j = <f>j + 1 — 0y- 1> 00 = —<*<j >2 + (1 + 2a)<£i — (1 + ot)(f> o. (5.5a,b) 

c c 

Since (5.5a) is a difference equation for we look for a solution of (5.5) of the form <f>j = /c ; and 
obtain 

= k , = (/c — 1 )[— ok + (1 + o)]. (5.6a, b) 

C K C 

By eliminating sAx/c between (5.6a, b), one obtains the following cubic equation for k: 

q(ii) = (k — l) 2 (2a/c — 1) = 0. (5.7) 

The roots of (5.7), i.e., the zeros of g(t c), are 


K — 1,1, 


2a 


(a f 0). 


(5.8) 


We then check to see if there is a GKS eigenvalue or a GKS generalized eigenvalue. We omit 
the details of the analysis (see [8]) and summarize the GKS stability results as follows: 


a < —1/2 unstable (GKS eigenvalue) 

a = -1/2 unstable (GKS generalized eigenvalue) 

a > — 1/2 stable (no GKS eigenvalue or generalized eigenvalue). 


(5.9a) 

(5.9b) 

(5.9c) 


6. Normal Mode Analysis (Finite-Domain Problem) 

The normal mode analysis for the finite domain proceeds in the same way as the right-quarter- 
plane normal mode analysis except that there are two boundaries and, consequently, the analytical 
boundary condition (3.2) is imposed at i = L. The results of the analysis (see [8]) are summarized 
as follows. An eigenfunction is 

^ = ^ + (-l) J+1 /c 2J (-l//c)^ (6.1) 

and the corresponding eigenvalue is 

1 A Ax /y% 

2,3 — K , S = -r-rS (6.2) 

K \C\ 


I 
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where k is a zero of the polynomial /(«) defined by 

f(ic) = (k — l) 2 (2a/c — 1) + k? j [- (— 1) J k~ 1 (k + 1) 2 (* + 2a)]. (6.3) 

A complete solution of the eigenvalue problem requires that we find the zeros of the polynomial 
(6.3). Unfortunately, one cannot obtain the zeros of /(k) analytically and so the eigenvalue problem 
for the finite domain is intractable. On the other hand, a knowledge of the precise values of the 
zeros is more information than we need since the algebraic tests for stability require only asymptotic 
(J — » oo) values of the eigenvalues. 

7. Finite Domain vs. Quarter Plane: the Eigenvalue Connection 
The /c-polynomial (6.3) for the finite domain problem can be written as 

/(«) = q(i c) + K? J h(K) (6.4) 


where 

g(*c) = (k — l) 2 (2a*c — 1), H(k) = [— (— l) J /c -1 (« + l) 2 (/c + 2a)]. (6.5a, b) 

We have intentionally split the polynomial /(/c) into two parts where the polynomial g(/c) is 
precisely the ^-polynomial (5.7) associated with the right-quarter-plane problem of the GKS theory. 

The notion of stability for the finite domain problem is intimately associated with the solution 
behavior as the spatial mesh is refined, i.e., J —* oo. Hence we are primarily interested in the zeros 
of /(/c) for large J. In particular, we are interested in the conditions under which the polynomial 
/(«) reduces to the quarter-plane polynomial q[n) in the limit J — ♦ oo. The reduction obviously 
occurs when |/c| 2,7 — ► 0 as J — ► oo. The resulting conditions (i.e., the asymptotic behavior of k) 
determine the connection between the normal mode analysis of the finite domain problem and the 
normal mode analysis of the quarter plane problem. 

If a zero of the polynomial /(«) can be estimated asymptotically (for large J), then the cor- 
responding eigenvalue s is given by (6.2). One can restate the necessary condition (4.2) for Lax- 
Richtmyer stability by requiring that every eigenvalue s of the matrix A satisfy 

*M < 7- (6-6) 

Lax-Richtmyer instability occurs if inequality (6.6) is not satisfied. In the asymptotic analysis to 
follow, we relate the presence of a GKS eigenvalue or GKS generalized eigenvalue to Lax-Richtmyer 
instability. 

What values of k do we actually need to consider? One can show that there is no loss in 
generality in assuming that |/c| < 1. In general, the zeros of f(n) which are crucial to the stability 
(actually instability) of a semidiscrete approximation depend on J and we write k = k(J). Since 
we are assuming |/c| < 1, we let 


|k| = |k(J)|=: 1- €) 0 < e(J) < 1. 

There are three possible cases to consider: 

(6.7) 

case (1) 

: e(J) > * >0 as J — ► oo 

(6.8a) 

case (2) 

: e(J) —* 0 as J — ♦ oo 

(6.8b) 

case (3) 

: e(J) = 0 for all J. 

(6.8c) 


From (6.2) it easy to show that the real part of s can be written as 

«(M) = *(« - l/«) = a(l "j) ~ 11 (6.9) 


5 


where k — a + ib, |«| 2 — a, 3 +b 3 . In case (l) SR(s) ^ 0 as J — ► oo, in case (2) !R(s) —► 0 as J — » oo, 
and in case (3) S?(s) = 0 for all J. 

case (1) 

In the first case we assume that e is strictly bounded away from zero as J — » oo, i.e. |/c| < 1 as 
J — ► oo and consequently 

lim |k| 2J = 0. (6.10) 

J— »oo 

Then in the limit J — ► oo, the polynomial /(«) reduces to the right-quarter-plane polynomial 

g(/c) = (« — l) a (2a/c — 1). (6.11) 

We next check to see if the above cubic polynomial has one or more zeros with |k| < 1. If there is 
no zero |k| < 1, then the assumption which led to (6.11) is invalid and we drop consideration of 
this case. A necessary condition for Lax-Richtmyer stability is that all the eigenvalues must satisfy 
inequality (6.6) and hence 8?(s) < 0 in the limit J — ► oo. But in obtaining (6.11) we have already 
taken the limit J — ► oo and the semidiscrete approximation will clearly be unstable if 9?(!) > 0. 
Consequently, if q(/c) = 0 and (6.2), i.e., the resolvent equations, have a nontrivial solution of the 
form 

iR(S) > 0, |/c| < 1, (6.12) 

then the semidiscrete approximation is Lax-Richtmyer unstable. It is clear that (6.12) is identical to 
(5.3) of the GKS theory. Consequently, case (1) with SR(S) > 0 corresponds to a GKS eigensolution. 

case (2) 

In case (1), |k| was assumed to be strictly less than unity as J — * oo. The second case of interest 
is |k| < 1 for any finite J but |/c| — ► 1~ as J — ► oo. From (6.7) and (6.8b) we have 

|k| = 1 — e(J), e(J) — * 0 as J — ► oo, (6.13) 


and consequently 


|«| 2J = [l-£(J)] 2J =e 3J,,, > 1 - <(J )l«e- 2J ‘( J ), oo. (6.14) 

Considering the product J • e(J) in the limit J — ► oo there are only two possibilities: either 

J ■ e(J) < K, J — ► oo, and lim |/c| 2J = constant > 0 (6.15a) 

where K is a positive constant independent of J, or 

J • e(J) — * oo, J — ► oo, and lim IkI 2-7 = 0. (6.15b) 

J~*oo 

In general, S?(s) is a function of k, e.g., (6.9), and using (6.13) we obtain 

5R(s) « We, J ~+ oo (6.16) 


where W is a constant. 

First we consider the possibility given by (6.15a). The /c-polynomial /(/c) does not reduce to the 
quarter-plane polynomial g(/c) , however 


3?(l) < 


\W\K 

J 


which satisfies the necessary condition (6.6) for Lax-Richtmyer stability. 
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For the second possibility given by (6.15b) the ^-polynomial /(/c) does, in fact, reduce to the 
quarter-plane polynomial qfa). But from (6.15b) and (6.16) we have for any finite value of K , and 
J sufficiently large, 

3t(s) > if W > 0 (6.17) 

which violates the necessary Lax-Richtmyer stability condition (6.6) . 

In summary, if the reduced k polynomial <?(/c) has a zero |/c| = 1 and a perturbation of the k 
inside the unit circle (|/c| = 1 — e, e > 0) leads to a positive S?(5) (W > 0 in (6.16)) the scheme is 
Lax-Richtmyer unstable. This algebraic test sequence is precisely the GKS test for a generalized 
eigenvalue (see (5.4))! 
case (3) 

This final case has the distinct feature that |«| = 1 is a solution of the k polynomial for all J and 
consequently I = 0 for all J. The k polynomial does not reduce to the quarter-plane polynomial. 
The GKS analysis does not distinguish this case from case (2) and both are considered in the GKS 
generalized eigenvalue test. From the point of view of an eigenvalue analysis, the two cases must 
be treated separately since any instability in the present case derives not from an eigenvalue with 
a positive real part but from the algebraic growth (as J — * oo) of the norm of the solution due 
to the normal mode with eigenvalue s = 0. Details of the analysis for this case will be presented 
elsewhere [8]. 

8. Concluding Remarks 

We have presented a brief outline of an analysis that correlates Lax-Richtmyer stability and the 
GKS normal mode stability analysis. A more complete analysis leads us to conjecture that, except 
for certain identifiable borderline cases 

GKS stability <=> Lax-Richtmyer stability, 
i.e., the GKS algebraic tests can be used to check Lax-Richtmyer stability. 
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